** Set filepaths to folders to store intermediate data files, output and .do files:
global datafiles "N:\SpecialProjects\Inequality_Egypt\Peru_2019\data"
global dofiles "N:\SpecialProjects\Inequality_Egypt\Peru_2019\do_files"

* The following .do file runs the estimation code for DHSV and populates a matrix with the estimates (before producing bootstrapped standard errors).  These are the estimates found in Table 4 (and B.3 and B.4).
do "$dofiles\DHIS_V_pre_boot_EDCC_replication_file.do"


use "$datafiles/peru_dhsV_analysis.dta", clear

**** The bootstrap routine starts here ****

/*** Set up a program called myboot to run the estimation for the bootstrap ***/
capture program drop myboot
	program define myboot, rclass
	preserve
	bsample
	  
	********************* Now execute the estimation procedure *********************************  
	
** Policy and analysis years for main specification (as seen in Table 4)
local aby "1990"
local aey "1997"
local pby "1996"
local pey "1997"
local exc "& year~=1995"

/*
** alternative specifications (Table 6):

** extend policy years from 1995 to 1998:
local aby "1990"
local aey "1998"
local pby "1995"
local pey "1998"
local exc " "

** extend policy years from 1995 to 1999:
local aby "1990"
local aey "1999"
local pby "1995"
local pey "1999"
local exc " "

*/

	gen bid=_n

	tempname bootsample
	save `bootsample', replace

************* Logit (as described in paper Section 6 and shown in Appendix Table B.1) ******************************

keep bid age_* ster_* kids_* mar_* will_end*  boys_*  born_* m_age_fb yr_st ed_cat weight1 rur_coast urb_mount rur_mount urb_jungle rur_jungle urb_coast 

*Reshape to create a pseudo-panel
reshape long age_ ster_ kids_  mar_ will_end_  boys_  born_  , i(bid) j(year)

* dummy for being sterilized in a given year
gen sterilized=yr_st==year


* Drop observations after sterilization (no longer at risk of sterilization)
bys bid: gen flag=cond(sterilized==0, sum(sterilized),0)
drop if flag==1


* It was "illegal" to be sterilized if you were not married and did not have kids - and almost noone was, so we make this 
* sample restriction (otherwise the probit does not converge due to perfect prediction of failure.)
gen eligable=(mar_ | will_end_) & kids_>0


gen kids_count=kids_
replace kids_count=. if kids_==0
replace kids_count=10 if kids_>10
tab kids_count, gen(kids)

gen year_lin=year-1992

gen young=age_<26
gen mid1_age=age_>=26 & age_<=29
gen mid2_age=age_>29 & age_<=36
gen older=age_>36

gen ed_prim=ed_cat==2 | ed_cat==1
gen ed_second=ed_cat==3
gen ed_high=ed_cat==4

gen no_birth=born_==0

gen no_boys=boys_==0


local int_vars " rur_mount rur_jungle ed_prim young mid1_age no_birth "

cap drop policy
gen policy=0
replace policy=1 if year>=`pby' & year<=`pey'


cap drop policy_X_*
foreach x of local int_vars{
qui gen policy_X_`x'=`x'*policy
}

forval i=1/9{
gen policy_X_kids`i'=policy*kids`i'
}


gen policy_save=.
foreach x of local int_vars{
qui gen policy_X_`x'_save=.
}

forval i=1/9{
gen policy_X_kids`i'_save=.
}


logit sterilized  kids1-kids9 m_age_fb rur_coast urb_mount rur_mount urb_jungle rur_jungle ed_prim ed_second young mid1_age older  no_boys no_birth  ///
      policy_X_rur_mount policy_X_rur_jungle policy_X_ed_prim policy_X_young policy_X_mid1_age  policy_X_no_birth  year_lin [pw=weight1] if year>=`aby' & year<=`aey' `exc' & eligable, vce(cluster bid)


qui cap drop pr2_

qui predict pr2_ if year>=`pby' & year<=`pey' `exc' & eligable & e(sample)


qui replace policy_save=policy
foreach x of local int_vars{
qui replace policy_X_`x'_save=policy_X_`x'
}

forval i=1/9{
qui replace policy_X_kids`i'_save=policy_X_kids`i'
}


qui replace policy=0
foreach x of local int_vars{
qui replace policy_X_`x'=0
}

forval i=1/9{
qui replace policy_X_kids`i'=0
}

qui cap drop pr1_


qui predict pr1_ if year>=`pby' & year<=`pey' `exc'  & eligable & e(sample)


qui replace policy=policy_save
foreach x of local int_vars{
qui replace policy_X_`x'=policy_X_`x'_save
}

forval i=1/9{
qui replace policy_X_kids`i'=policy_X_kids`i'_save
}

save "$datafiles/peru_dhsV_qp.dta", replace

** All de jure women who were eiligable in 1996 or 1997 have a pr1 pr2 for that year 
** Being eligible means having ever been married and having at least one child	

keep pr1_* pr2_*  bid year 

	* make this wide again
	reshape wide   pr1_ pr2_  , i(bid) j(year)
	sort bid
	save  "$datafiles/pscoresV.dta", replace

	** Merge the p1 p2 and delta_p variables onto the cross-sectional data set:
	use `bootsample', clear
	cap drop _merge
	sort bid
	merge bid using "$datafiles/pscoresV.dta"

cap drop S
gen S = yr_st>=`pby' & yr_st<=`pey' 


gen status=.
replace status=1 if yr_st<`pby' & yr_st~=.               /* sterilized before the policy -- these women are all excluded from the outcome results below */


egen _p2=rowmean(pr2_*)
egen _p1=rowmean(pr1_*)

gen delta_p=_p2-_p1
gen dp_p2=delta_p/_p2

*** Standard Reweighting 

gen sdrweight=weight1 if S==1 & status~=1
replace sdrweight=weight1*(_p2/(1-_p2)) if S==0 & status~=1

bys S: egen sdsrweight_t=total(sdrweight)
gen sdsrweight=sdrweight/sdsrweight_t


*** Our reweighting S2=1, C=1
gen rweight=weight1*dp_p2 if S==1 & status~=1                 
replace rweight=weight1*(delta_p/(1-_p2)) if S==0 & status~=1 

bys S: egen srweight_t=total(rweight)
gen srweight=rweight/srweight_t

*** Our reweighting S2=1, C=0

gen orweight=weight1*(_p1/_p2) if S==1 & status~=1           
replace orweight=weight1*((_p1)/(1-_p2)) if S==0 & status~=1 

bys S: egen osrweight_t=total(orweight)
gen osrweight=orweight/osrweight_t



	tempname bootsample_ww
	save `bootsample_ww', replace
	
**************************** Estimating Results *********************************************************************************
	
	******** estimating the number of women affected *********************
	gen w3=weight1*325.8
	total w3 if S==1
	matrix total=e(b)
	return scalar women_all=total[1,1]

	gen estim=w3*dp_p2
	total estim if S==1
	return scalar women_campaign=_b[estim]


	**************************** Number of Kids *****************************************************************
	* Modified Reweighting (S2=1, C=1) 
	reg kids_2004 S  [pw=srweight] if status~=1 
	return scalar kids_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg kids_2004 S [pw=osrweight] if status~=1 
	return scalar kids_s2d0=_b[S]



	**************************** Mother's Labor-Force Participation (working for pay) ***************************


	* Modified Reweighting (S2=1, C=1) 
	reg working_p S  [pw=srweight] if status~=1 
	return scalar wp_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg working_p S [pw=osrweight] if status~=1 
	return scalar wp_s2d0=_b[S]


	**************************** Mother's Health *****************************************************************
				 ** weight_height_m

	* Modified Reweighting (S2=1, C=1) 
	reg weight_height_m S  [pw=srweight] if status~=1 & pregnant==0
	return scalar WH_M_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg weight_height_m S [pw=osrweight] if status~=1  & pregnant==0
	return scalar WH_M_s2d0=_b[S]

				 ** height_age_m

	* Modified Reweighting (S2=1, C=1) 
	reg height_age_m S  [pw=srweight] if status~=1 & pregnant==0
	return scalar HA_M_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	reg height_age_m S [pw=osrweight] if status~=1  & pregnant==0
	return scalar HA_M_s2d0=_b[S]

				   **Anemic

	* Modified Reweighting (S2=1, C=1) 
	reg anemic_m S  [pw=srweight] if status~=1 & pregnant==0
	return scalar anemic_M_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg anemic_m S [pw=osrweight] if status~=1 & pregnant==0
	return scalar anemic_M_s2d0=_b[S]


	*********** Domestic Violence LAST 12 MONTHS*******************************************************
	
	* Modified Reweighting (S2=1, C=1) 
	reg p_physexviol12m S [pw=srweight] if status~=1
	return scalar dm_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	reg p_physexviol12m S [pw=osrweight] if status~=1 
	return scalar dm_s2d0=_b[S]


	***********EDUCATION OF CHILDREN

	**********Education of household children (under 15) but born prior to the program
	** need to use the household record to get these kids' education which leads to slight differences

	use `bootsample_ww', clear

	
	keep kid_age_ed* kbyed* kid_boy_ed* enrolled* sch_yr_* grade_at_age* bid status S srweight osrweight

	reshape long kid_age_ed kbyed kid_boy_ed enrolled sch_yr_ grade_at_age, i(bid) j(kid)

	rename kid_boy_ed kid_boy
	rename kid_age_ed kid_age
	rename kbyed kby




	*********** Education of own children**********************************************************************


	*********** Enrollment of own kids**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 
	return scalar ENkid_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg enrolled S i.kid_age [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 
	return scalar ENkid_s2d0=_b[S]

	gen ed_flag=e(sample)


	*********** grade for age of own children**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6 & ed_flag==1
	return scalar GFAkid_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & ed_flag==1
	return scalar GFAkid_s2d0=_b[S]



	**********************ED by GENDER ***************************************************************************


	*********** Enrollment of own boys**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==1
	return scalar ENboy_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg enrolled S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==1 
	return scalar ENboy_s2d0=_b[S]

	*********** Enrollment of own girls**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==0
	return scalar ENgirl_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg enrolled S i.kid_age [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==0 
	return scalar ENgirl_s2d0=_b[S]



	********** grade for age of own boys **********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==1
	return scalar GFAboy_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==1
	return scalar GFAboy_s2d0=_b[S]

	********** grade for age of own girls **********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==0
	return scalar GFAgirl_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==0
	return scalar GFAgirl_s2d0=_b[S]




	restore
end



simulate    kids1=r(kids_s2d1)	kids0=r(kids_s2d0) WHM1=r(WH_M_s2d1) WHM0=r(WH_M_s2d0) HAM1=r(HA_M_s2d1) HAM0=r(HA_M_s2d0) ancM1=r(anemic_M_s2d1) ancM0=r(anemic_M_s2d0)  DM1=r(dm_s2d1) DM0=r(dm_s2d0) wp1=r(wp_s2d1)	wp0=r(wp_s2d0)	ENkid1=r(ENkid_s2d1) ENkid0=r(ENkid_s2d0) 	GFA1=r(GFAkid_s2d1) GFA0=r(GFAkid_s2d0)	 ENboy1=r(ENboy_s2d1)	ENboy0=r(ENboy_s2d0) ENgirl1=r(ENgirl_s2d1)	ENgirl0=r(ENgirl_s2d0) GFAboy1=r(GFAboy_s2d1) GFAboy0=r(GFAboy_s2d0) GFAgirl1=r(GFAgirl_s2d1) GFAgirl0=r(GFAgirl_s2d0)	women_all=r(women_all)	 women_campaign=r(women_campaign) , reps(500) seed(8675309): myboot

bstat, stat(non_boot)
